Feature Angle and Dihedral Restraints #272

akalpokas · 2025-01-10T17:56:24Z

Introduction

This PR expands the restraints functionality of sire in order to include stand-alone angle and dihedral restraints, which can be used in normal and alchemical simulations within sire. This restraints expansion in theory enables scaffold-hopping & ring-breaking FEP within sire through the use of auxiliary restraints approach (https://pubs.acs.org/doi/10.1021/acs.jctc.1c00214), although this has only been tested for cyclohexane --> hexane transformation. Regardless of whether this is the actual approach we will proceed in the future for scaffold-hopping experiments, expansion of restraints functionality will allow for cross-comparison with other scaffold-hopping methodologies, and further expands the capabilities of sire dynamics.

Proposed Changes

Specifically this PR does the following:

Adds support for angle and dihedral restraints in the C++ layer (as part of SireMM corelib). Note that this replaces previously unused boilerplate angle and dihedral restraint files located there. Most of this code is adapted and modified from already existing restraints classes.
Updates sire to openMM system code (sire_to_openmm_system.cpp) to create a CustomAngleForce and a CustomTorsionForce based on these restraints. Note that currently the angle and torsion restraint forces are specified and applied in a harmonic functional form. This was done in order to keep consistent with Boresch restraint force implementation, but can be easily changed if needed.
Updates Lambda lever code (lambdalever.cpp) to allow scaling of the restraints force with rho - this effectively allows these restraints to be used in alchemical simulations just as any other perturbable restraint.
Builds python wrappers for all of the modifications above.
Exposes these restraints to the sire python API layer.
Expands the restraints tutorial showcasing the creation of these restraints for an alanine molecule, and adds basic unit testing for the python code.

In the end the user can create angle and dihedral restraints in the same manner as other already existing restraints. For example:

import sire as sr
mols = sr.load(sr.expand(sr.tutorial_url, "ala.top", "ala.crd"))
restraints = sr.restraints.angle(mols=mols, atoms=[0, 1, 2])
print(restraints)
AngleRestraints( name=restraint, size=1
0: AngleRestraint( [0, 1, 2], theta0=112.8°, ktheta=0.0304617 kcal mol-1 °-2 )
 )

or

import sire as sr
mols = sr.load(sr.expand(sr.tutorial_url, "ala.top", "ala.crd"))
restraints = sr.restraints.dihedral(mols=mols, atoms=[0, 1, 2, 3])
print(restraints)
DihedralRestraints( name=restraint, size=1
0: DihedralRestraint( [0, 1, 2, 3], phi0=244.528°, kphi=0.0304617 kcal mol-1 °-2 )
)

Because most of the restraints code is adapted from already existing restraints classes, these restraints support the same functionality as others (custom force values, usage of containers for passing atoms, multiple restraints, etc.)

Example Usage

This section showcases angle and dihedral restraints in action for alchemical simulations for a hexane molecule. Please compile the code from this PR as well as BioSimSpace in order to run this code.

Alchemical Angle Restraints

In this example, we are going to restrain the first angle in a hexane molecule to 3 rad by applying a force of 100 kcal mol-1 rad-2 by switching restraints force as a function of lambda (so that restraints are fully on at lambda=1). Please create alchemical_angle_restraints_test directory before running this test.

Setup and Dynamics

import sire as sr
import BioSimSpace as BSS

import mdtraj as md
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Build a hexane molecule
hexane = BSS.Parameters.gaff2("CCCCCC").getMolecule()
mols = BSS.Convert.toSire(hexane)

# Create a lambda schedule where the restraint lever is linearly scaled from 0 to 1 with lambda
l = sr.cas.LambdaSchedule()
l.add_stage("morph", l.initial())
l.set_equation(stage="morph", lever="restraint", equation=l.initial() * l.lam())

# Define the angle restraints
ang = mols.angles()[0]
restraints = sr.restraints.angle(mols=mols, atoms=ang.atoms(), theta0="3 rad", ktheta="100 kcal mol-1 rad-2")

# Setup the dynamics
lambda_values = [x / 100.0 for x in range(0, 101, 10)]

timestep = "2fs"
energy_frequency = "100ps"
constraint = "h-bonds"
perturbable_constraint = None
equil_time = "5ps"
run_time = "20ps"

mols = mols.minimisation().run().commit()

# Run the dynamics
for i, lambda_value in enumerate(lambda_values):
    print(f"Simulating lambda={lambda_value:.2f}")
    min_mols = (
        mols.minimisation(
            lambda_value=lambda_value,
            constraint=constraint,
            perturbable_constraint="none",
        )
        .run()
        .commit()
    )

    d = mols.dynamics(timestep="2fs",
                            temperature="25oC",
                            restraints=restraints,
                            schedule=l,
                            constraint=constraint,
                            perturbable_constraint=perturbable_constraint,
                            lambda_value=lambda_value,
                            )

    d.randomise_velocities()
    d.run(equil_time, save_frequency=0)
    print("Equilibration complete")

    d.run(
        run_time,
        energy_frequency=energy_frequency,
        frame_frequency="1ps",
        lambda_windows=lambda_values,
    )
    print("Dynamics complete")
    d.to_xml(f"alchemical_angle_restraints_test/debug_alchemical_restraints_{i}.xml")
    sr.save(d.commit().trajectory(), f"alchemical_angle_restraints_test/alchemical_restraints_test_{i}", format=["PRM7", "DCD"])
    mols.delete_all_frames()

Plotting

After the dynamics have concluded, we can visualize the restraint being turned on:

df_restraints_combined = []

for i in range(0,11):
    traj = md.load(f"alchemical_angle_restraints_test/alchemical_restraints_test_{i}.dcd", top=f"alchemical_angle_restraints_test/alchemical_restraints_test_{i}.prm7")
    top = traj.topology
    c1 = top.select("name C1")
    c2 = top.select("name C2")
    c3 = top.select("name C3")
    c4 = top.select("name C4")

    sel0 = np.array([c1, c2, c3])
    sel0 = sel0.reshape(-1, 3)
    angles = md.compute_angles(traj, sel0)
    frame = traj.time


    df_restraints = pd.DataFrame()
    df_restraints["angles"] = angles.flatten()

    df_restraints["frame"] = frame
    df_restraints["lambda_window"] = i
    df_restraints_combined.append(df_restraints)

df_restraints_combined = pd.concat(df_restraints_combined)

fig, ax = plt.subplots(figsize=(8,8))
sns.lineplot(data=df_restraints_combined, x="frame", y="angles", hue="lambda_window")
ax.set_title("Angle Between C1 C2 and C3")
ax.set_ylabel("Angle (radians)")
ax.set_xlabel("Frame")
plt.savefig("angle_restraints.png")

Which gives:

We can see that as lambda is increasing, the average value of the first angle is increasing from around 1.9 rad (equilibrium value is 1.96 rad for this particular angle) to around 2.55 rad.

We can also inspect the openMM system files created in alchemical_angle_restraints_test directory to see the CustomAngleForce. E.g. for the first lambda window:

		</Force>
		<Force energy="rho*k*(theta-theta0)^2;" forceGroup="0" name="CustomAngleForce" type="CustomAngleForce" usesPeriodic="1" version="3">
			<PerAngleParameters>
				<Parameter name="rho"/>
				<Parameter name="k"/>
				<Parameter name="theta0"/>
			</PerAngleParameters>
			<GlobalParameters/>
			<EnergyParameterDerivatives/>
			<Angles>
				<Angle p1="0" p2="1" p3="2" param1="0" param2="418.40000000000003" param3="3"/>
			</Angles>
		</Force>

Alchemical Dihedral Restraints

In this example, we are going to restrain the first torsion in a hexane molecule to 2 rad by applying a force of 10 kcal mol-1 rad-2 by again switching restraints force as a function of lambda. Please create alchemical_dihedral_restraints_test directory before running this test.

Setup and Dynamics

import sire as sr
import BioSimSpace as BSS

import mdtraj as md
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Build a hexane molecule
hexane = BSS.Parameters.gaff2("CCCCCC").getMolecule()
mols = BSS.Convert.toSire(hexane)

# Create a lambda schedule where the restraint lever is linearly scaled from 0 to 1 with lambda
l = sr.cas.LambdaSchedule()
l.add_stage("morph", l.initial())
l.set_equation(stage="morph", lever="restraint", equation=l.initial() * l.lam())

# Define the angle restraints
dih = mols.dihedrals()[0]
restraints = sr.restraints.dihedral(mols=mols, atoms=dih.atoms(), phi0="2 rad", kphi="10 kcal mol-1 rad-2")

# Setup the dynamics
lambda_values = [x / 100.0 for x in range(0, 101, 10)]

timestep = "2fs"
energy_frequency = "100ps"
constraint = "h-bonds"
perturbable_constraint = None
equil_time = "5ps"
run_time = "20ps"

mols = mols.minimisation().run().commit()

# Run the dynamics
for i, lambda_value in enumerate(lambda_values):
    print(f"Simulating lambda={lambda_value:.2f}")
    min_mols = (
        mols.minimisation(
            lambda_value=lambda_value,
            constraint=constraint,
            perturbable_constraint="none",
        )
        .run()
        .commit()
    )

    d = mols.dynamics(timestep="2fs",
                            temperature="25oC",
                            restraints=restraints,
                            schedule=l,
                            constraint=constraint,
                            perturbable_constraint=perturbable_constraint,
                            lambda_value=lambda_value,
                            )

    d.randomise_velocities()
    d.run(equil_time, save_frequency=0)
    print("Equilibration complete")

    d.run(
        run_time,
        energy_frequency=energy_frequency,
        frame_frequency="1ps",
        lambda_windows=lambda_values,
    )
    print("Dynamics complete")
    d.to_xml(f"alchemical_dihedral_restraints_test/debug_alchemical_restraints_{i}.xml")
    sr.save(d.commit().trajectory(), f"alchemical_dihedral_restraints_test/alchemical_restraints_test_{i}", format=["PRM7", "DCD"])
    mols.delete_all_frames()

Plotting

df_restraints_combined = []

for i in range(0,11):
    traj = md.load(f"alchemical_dihedral_restraints_test/alchemical_restraints_test_{i}.dcd", top=f"alchemical_dihedral_restraints_test/alchemical_restraints_test_{i}.prm7")
    top = traj.topology
    c1 = top.select("name C1")
    c2 = top.select("name C2")
    c3 = top.select("name C3")
    c4 = top.select("name C4")

    sel0 = np.array([c1, c2, c3, c4])
    sel0 = sel0.reshape(-1, 4)
    dihedrals = md.compute_dihedrals(traj, sel0)
    frame = traj.time


    df_restraints = pd.DataFrame()
    df_restraints["dihedrals"] = dihedrals.flatten()

    df_restraints["frame"] = frame
    df_restraints["lambda_window"] = i
    df_restraints_combined.append(df_restraints)

df_restraints_combined = pd.concat(df_restraints_combined)

fig, ax = plt.subplots(figsize=(8,8))
sns.lineplot(data=df_restraints_combined, x="frame", y="dihedrals", hue="lambda_window")
ax.set_title("Dihedral Between C1 C2 C3 and C4")
ax.set_ylabel("dihedral (radians)")
ax.set_xlabel("Frame")
plt.savefig("dihedral_restraints.png")

Which gives:

We can see that as the value of lambda increases, the dihedral gets restrained to a value of around 2 rad.

Again we can inspect the .xml file to see the force in action:

		</Force>
		<Force energy="rho*k*min(dtheta, two_pi-dtheta)^2;dtheta = abs(theta-theta0);two_pi=6.283185307179586;" forceGroup="0" name="CustomTorsionForce" type="CustomTorsionForce" usesPeriodic="1" version="3">
			<PerTorsionParameters>
				<Parameter name="rho"/>
				<Parameter name="k"/>
				<Parameter name="theta0"/>
			</PerTorsionParameters>
			<GlobalParameters/>
			<EnergyParameterDerivatives/>
			<Torsions>
				<Torsion p1="0" p2="1" p3="2" p4="3" param1="0" param2="41.84" param3="2"/>
			</Torsions>
		</Force>

Checklists

I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): [y]
I confirm that I have added a test for any new functionality in this pull request: [y]
I confirm that I have added documentation (e.g. a new tutorial page or detailed guide) for any new functionality in this pull request: [y]
I confirm that I have added a changelog entry to the changelog (we will add a link to this PR as part of the review): [y]
I confirm that I have permission to release this code under the GPL3 license: [y]

Suggested reviewers:

@lohedges

Happy to make further improvements if needed!

…s update

…nMM system

… bugs might still be lurking somewhere

… also renamed the core alchemical angle/dihedral restraint c++ files

… also slightly tidied up python API for these restraints

lohedges · 2025-01-13T09:18:35Z

Thanks, @akalpokas. This looks great! I'll go through this during the week in between doing other things. I'll probably do it in chunks, so expect to see comments/suggestions drop in now and then.

src/sire/restraints/_restraints.py

wrapper/Convert/SireOpenMM/lambdalever.cpp

wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp

…es [ci skip]

lohedges · 2025-01-17T15:56:49Z

Thanks, @akalpokas, this all looks good. My only thoughts are regarding the functional forms and if we want to support anything different in future. For example, we've had some recent interest in the restricted bending potential from GROMACS. I'm not sure if it makes sense to leave harmonic potentials as the default for angle and dihedral, then give more exotic ones other names, or support something like angle.harmonic, or rename to harmonicangle, etc. I'll have a thought over the weekend and let you know.

akalpokas · 2025-01-17T16:27:02Z

Hi @lohedges, thanks for your review. I agree that the harmonic potential might not be ideal long term solution. When I was implementing this originally I was wondering if I could implement the restraint in such a way that the functional form would be flexible, i.e. the user could specify the expression for the OpenMM force upon creation of the restraint as one of the arguments. I think it's possible to modify the restraints class so that it has OpenMM expression as a parameter upon construction, sets it to harmonic functional form as default or takes the user specific one instead and uses it during Sire to OpenMM system creation.

corelib/src/libs/SireMM/dihedralrestraints.cpp

lohedges · 2025-01-17T16:52:14Z

corelib/src/libs/SireMM/dihedralrestraints.cpp

+    // if (distinct.count() != 3)
+    //     throw SireError::invalid_arg(QObject::tr(
+    //                                      "There is something wrong with the atoms provided. "
+    //                                      "They should all be unique and all greater than or equal to 0."),
+    //                                  CODELOC);


More commented code.

corelib/src/libs/SireMM/anglerestraints.cpp

lohedges · 2025-01-17T16:54:14Z

I think that could work, but might be in tricky in general since we'd need a way of coupling a parameter in the expression to what's used by the lambda lever. Will have a think.

akalpokas · 2025-01-21T15:53:17Z

@lohedges I have addressed the suggested changes, thanks for the review. Happy for this to be merged!

lohedges

Many thanks for this, I'll merge across to my branch and merge to devel in one-shot.

corelib/src/libs/SireMM/anglerestraints.cpp

akalpokas added 14 commits July 10, 2024 11:14

Initial update for angle restraints

05e38a9

Intial addition of alchemical dihedral restraints and angle restraint…

12af94a

…s update

Update wrappers for angle and dihedral restraints classes

8809ac5

Added intial python code for implementing alchemical dihedral restraints

386912d

Added prototype code for adding alchemical dihedral restraints to Ope…

69a7ee2

…nMM system

Fix conversion of dihedral restraints to OpenMM system

c29c4eb

Added full end-to-end alchemical dihedral restraints implementation -…

4a0726d

… bugs might still be lurking somewhere

Added full end-to-end alchemical angle restraints implementation, and…

9eeb824

… also renamed the core alchemical angle/dihedral restraint c++ files

Add the missing corelib files for angle restraints

a500bc9

Merge branch 'OpenBioSim:devel' into feature_alchemical_restraints_contd

55eb13e

Sync with upstream devel branch

c5404cf

Updated restraints tutorial to include angle and dihedral restraints,…

691223a

… also slightly tidied up python API for these restraints

Add basic unit tests for the new restraints

ab2188d

Updated changelog and contributor pages for the PR

35725d8

akalpokas temporarily deployed to sire-build January 10, 2025 18:12 — with GitHub Actions Inactive

akalpokas had a problem deploying to sire-build January 10, 2025 18:12 — with GitHub Actions Failure